Bootstrap Inferenece

Edward Vytlacil

Motivation: Mincer Wage Equation

  • Motivation: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Implementing Bootstrap: Average Effect of Experience

  • Implementing Bootstrap: Experience level that maximizes log wages.

  • Other Bootstrap Variants

Motivation: Mincer Wage Equation

Mincer wage equation (Mincer 1958):

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

  • For ease of exposition, dropping education.

Motivation: Mincer Wage Equation

Mincer wage equation (Mincer 1958):

\[\ln(wage)_i = \beta_0 + \beta_1 \mbox{exp}_i + \beta_2 \mbox{exp}_i^2 + e_i\]

  • Average effect of experience: \[\beta_1 + 2 \cdot \beta_2 \cdot \mathbb{E}(\mbox{exp}),\]

  • Experience level that maximizes log wage: \[\left | \frac{\beta_1}{2 \beta_2}\right |.\]

Motivation: Mincer Wage Equation

  • Let \(\widehat \beta_{1,N}, \widehat \beta_{2,N}\) denote OLS estimators,

  • Let \(\overline{\mbox{exp}}_N\) denote sample mean of experience.

  • Estimate average effect of experience by \[\widehat \beta_{1,N} + 2 \cdot \widehat \beta_{2,N} \cdot \overline{\mbox{exp}}_N,\]

    • Estimator is a non-linear function of \((\widehat \beta_{1,N}, \widehat \beta_{2,N}, \overline{\mbox{exp}}_N)\).

Motivation: Mincer Wage Equation

  • Let \(\widehat \beta_{1,N}, \widehat \beta_{2,N}\) denote OLS estimators,

  • Estimate experience level that maximizes log wage by \[\left | \frac{\widehat \beta_{1,N}}{2 \widehat \beta_{2,N}}\right|.\]

    • Estimator is a non-linear functions of \((\widehat \beta_{1,N}, \widehat \beta_{2,N})\).

Implement Mincer Wage Equation

library (AER)
data(CPS1985)
df  <- CPS1985[ CPS1985$education==12, c("wage","experience")]
mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
  • Consider estimation using CPS 1985 data.

    • available in package AER.

    • after loading AER, see ?CPS1985 for description of data.

  • Consider only using observations with \(12\) years of schooling.

Implement Mincer Wage Equation

library (AER)
data(CPS1985)
df  <- CPS1985[ CPS1985$education==12, c("wage","experience")]
mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
[1] "Mean Experience: 18.14"

Call:
lm(formula = I(log(wage)) ~ experience + I(experience^2), data = df)

Coefficients:
    (Intercept)       experience  I(experience^2)  
      1.5829802        0.0352680       -0.0005831  

Implement Mincer Wage Equation

mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
avg.eff.exp
maxexp
[1] "Estimated Avg Effect Experience: 0.0141"
[1] "Estimated Experience that Maximizes Wage: 30.2"

Implement Mincer Wage Equation

mean.exp <-mean(df$experience) 
mean.exp
reg.1 <- lm(I(log(wage))~experience+I(experience^2),data=df)
reg.1
a <- c(0,1,2*mean.exp)
avg.eff.exp <- as.numeric(t(a)%*%coef(reg.1)) 
maxexp <- abs(coef(reg.1)[2]/(2*coef(reg.1)[3]))
avg.eff.exp
maxexp
  • \(\overline{\mbox{exp}}_N=\) 18.1

  • \(\widehat \beta_{1,N} =\) 0.0353

  • \(\widehat \beta_{2,N} =\) -0.000583

  • \(\widehat \beta_{1,N} + 2 \widehat \beta_{2,N} \overline{\mbox{exp}}_N=\) 0.0141

  • \(\left | \frac{\widehat \beta_{1,N}}{2 \widehat \beta_{2,N}}\right| =\) 30.2

Motivation: Mincer Wage Equation

  • \(\widehat \beta_{1,N} + 2 \widehat \beta_{2,N} \overline{\mbox{exp}}_N=\) 0.0141

  • \(\left | \frac{\widehat \beta_{1,N}}{2 \widehat \beta_{2,N}}\right| =\) 30.2

  • We have not covered methods for constructing standard errors or performing inference in this context.

  • How to compute standard errors? construct confidence intervals? perform hypothesis tests?

Delta Method vs Bootstrap

  • Motivation: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Implementing Bootstrap: Average Effect of Experience

  • Implementing Bootstrap: Experience level that maximizes log wages.

  • Other Bootstrap Variants

Traditional approach: Delta Method.

  • Asymptotic approximation based on linear approximation from Mean Value Theorem combined with LLN, CLT.

  • Limitations:

    • only justified asymptotically,
    • sometimes provides poor approximation,
    • derivations using delta method can be tedious, chance of algebraic errors.
  • we will not cover in Econ 123

Use Monte Carlo Simulation?

  • If knew distribution of population being sampled, could use Monte Carlo simulation drawing from that distribution to find distribution of our estimator.

    • But we don’t know population distribution.

Bootstrap principle

  • Treat sample of size \(N\) as if it were the population.

  • Perform MC simulation, creating many bootstrap samples of size \(N\) by resampling with replacement from original sample.

  • Compute estimator on each sample, use distribution of estimates across samples to approximate distribution of estimator.

  • Difference between drawing from population vs from sample should vanish as \(N \rightarrow \infty\).

Bootstrap: Basic Principle

Bootstrap Method

  • Limitations:

    • only justified asymptotically,
    • sometimes provide poor approximation,
    • can be computationally intensive
  • Advantages:

    • Can bypass derivations, complicated formulas. Less tedious to implement. Removes chance of algebraic error in derivations.
    • Sometimes provides better approximation than delta method.

Bootstrap Method

  • There are many versions of the bootstrap, but we will focus on nonparametric bootstrap and corresponding bootstrap:

    • estimate of bias,
    • standard errors,
    • normal-approximation CI,
    • percentile CI,
    • BCa: bootstrap accelerated bias-corrected percentile CI.

Bootstrap Algorithm

  • Motivation: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Implementing Bootstrap: Average Effect of Experience

  • Implementing Bootstrap: Experience level that maximizes log wages.

  • Other Bootstrap Variants

Sample and Estimator

  • Suppose \(\mathbf{W} = (W_1,W_2,...,W_N)\) is observed sample
    • e.g., \(W_i = (Y_i,X_i)\) in regression context.
  • Estimator \(\widehat \theta = s(\mathbf{W})\)
    • expressing estimator as function of data, rule for how to use data to compute estimate.
    • in our application, \(\widehat \theta\) is estimator of average effect of experience or of experience level that maximizes wages.

Bootstrap algorithm

Bootstrap Samples

Bootstrap Replications

  • Create \(B\) bootstrap samples, each constructed by drawing from original sample with replacement \(N\) times.

  • Compute bootstrap replication on each bootstrap sample, \(\widehat \theta^*(b) = s(\mathbf{W}^{*b})\)

\[\begin{align} \mathbf{W}^{*1} & = (W^{*1}_1,W^{*1}_2,…,W^{*1}_N) ~\Rightarrow \widehat \theta^*(1)=s(\mathbf{W}^{*1}) \\ \mathbf{W}^{*2} & = (W^{*2}_1,W^{*2}_2,…,W^{*2}_N) ~ \Rightarrow \widehat \theta^*(2)=s(\mathbf{W}^{*2}) \\ & \qquad \qquad \quad ~~~\vdots \qquad \qquad \qquad ~~~ \vdots \\ \mathbf{W}^{*B} & = (W^{*B}_1,W^{*B}_2,…,W^{*B}_N) \Rightarrow \widehat \theta^*(B)=s(\mathbf{W}^{*B}). \\ \end{align}\]
  • Use empirical distribution of \((\widehat \theta^*(1), \widehat \theta^*(2),...,\widehat\theta^*(B))\) to approximate distribution of \(\widehat \theta\).

Example: Algorithm for Bootstrap s.e.

  • Create \(B\) bootstrap samples, each constructed by drawing from original sample with replacement \(N\) times.

  • Compute bootstrap replication \(\widehat \theta^*(b)\) on each bootstrap sample.

  • Define bootstrap s.e. for \(\widehat \theta\) by \[\widehat{\mbox{se}}_B=\mbox{s.d.}(\widehat \theta^*(1),\widehat \theta^*(2),…,\widehat\theta^*(B))\] (illustrated in figure to right).

Normal Approx. Bootstrap CI

Given bootstrap standard error, \(\widehat{\mbox{se}}_B\):

  • Construct 95% normal-approximation bootstrap CI on \(\theta\) by \[C^{nb} = [\widehat{\theta}-1.96 \cdot \widehat{\mbox{se}}_B, ~\widehat{\theta}+1.96 \cdot \widehat{\mbox{se}}_B]\]

  • More generally, construct \(1-\alpha\) level normal-approximation bootstrap CI on \(\theta\) by \[[\widehat{\theta}- z_{1-\alpha/2} \cdot \widehat{\mbox{se}}_B, ~\widehat{\theta}+z_{1-\alpha/2} \cdot \widehat{\mbox{se}}_B]\] where \(z_{1-\alpha/2}\) is \(1-\alpha/2\) quantile of standard normal distribution.

Example: Bootstrap Estimate Bias

  • Given \(B\) bootstrap replications, define the boostrap estimate of bias by \[\widehat{\mbox{bias}}^*= \frac{1}{B} \sum_{b=1}^B \widehat \theta^*(b) - \widehat \theta \] where \(\widehat \theta\) is estimate from original sample.

  • why centered at \(\widehat \theta\)?

Example: Algorithm for Percentile CI

  • Given \(B\) bootstrap replications, define \(q^*_{\alpha}\) as the \(\alpha\) empirical quantile of \((\widehat \theta^*(1),\widehat \theta^*(2),…,\widehat\theta^*(B))\)

    • compute \(q^*_{\alpha}\) as the value such that fraction \(\alpha\) of the bootstrap replications are less than \(q^*_{\alpha}\).
  • Bootstrap Percentile CI: \[C^{PC}= [q_{\alpha/2},~ q_{1-\alpha/2}].\]

Pros and Cons Bootstrap Methods

  • Bootstrap standard errors:
    • Pros:
      • intuitive, easy to implement,
      • researchers familiar with s.e., is the norm to report s.e..
    • Cons:
      • Valid only under relatively strong conditions,
      • May be erratic, poorly behaved.

Pros and Cons Bootstrap Methods

  • Normal-approximation bootstrap CI:

    • Pros:
      • Intuitive, easy to implement,
      • Researchers familiar with normal CIs,
      • Many estimators are asymptotically normal.

Pros and Cons Bootstrap Methods

  • Normal-approximation bootstrap CI:

    • Cons:
      • Invalid if bootstrap s.e. are invalid,
      • Poorly behaved if bootstrap s.e. are poorly behaved,
      • Normality may be bad approximation
        • especially poor approximation if estimator is biased and/or has highly asymmetric distribution.

Pros and Cons Bootstrap Methods

  • bootstrap percentile CI:

    • Pros:

      • intuitive, easy to implement,

      • valid under weaker conditions than bootstrap s.e. or bootstrap normal CI.

    • Cons:

      • can provide poor approximation if estimator is biased and/or has highly asymmetric distribution.

Pros and Cons Bootstrap Methods

  • BCa:

    • “bootstrap accelerated bias-corrected percentile interval”

    • adjusts percentile CI for median bias and skewness in bootstrap distribution

Pros and Cons Bootstrap Methods

  • BCa:

    • Pros:
      • provides more accurate CIs than above options, especially when estimator is biased or has skewed distribution.
    • Cons:
      • complicated, not intuitive.
      • (somewhat) more computationally intensive.

Implementing Bootstrap: Average Effect of Experience

  • Motivation: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Implementing Bootstrap: Average Effect of Experience

  • Implementing Bootstrap: Experience level that maximizes log wages.

  • Other Bootstrap Variants

Implement: Avg Effect of Experience

  • We now implement bootstrap s.e. and CIs for average effect of experience.

    • \(\theta = \beta_{1} + 2 \cdot \beta_{2} \cdot \mathbb{E}[\mbox{exp}],\)

    • \(\widehat \theta = \widehat \beta_{1,N} + 2 \cdot \widehat \beta_{2,N} \cdot \overline{\mbox{exp}}_N.\)

  • We first proceed in base R, then consider using package boot.

  • First step: create bootstrap samples.

Implement: Avg Effect of Experience

  • Create bootstrap sample:

    • Randomly draw \(N\) integers, with replacement, from \((1,...,N)\), with each integer equally likely to be draw,

    • Create new sample formed by observations with those indices.

    • Some observations from original sample will be drawn multiple times into bootstrap sample, others will not appear in bootstrap sample.

N <- nrow(df) 
indices.b1<-sample(1:N,replace=TRUE)  # draw indices with replacement from (1,...,N), N sample size
tab.b1 <- df[indices.b1,] # construct bootstrap sample

Bootstrap Sample 1

Original Sample

Bootstrap Sample 1

Bootstrap Sample 2

Original Sample

Bootstrap Sample 2

Bootstrap Sample 3

Original Sample

Bootstrap Sample 3

Bootstrap Sample 4

Original Sample

Bootstrap Sample 4

Bootstrap Replications Avg. Effect Experience

  • given bootstrap sample, \(\mathbf{W}^{*b}=(W^{*b}_1,W^{*b}_2,…,W^{*b}_N)\), compute \(\widehat \theta^*(b)=s(\mathbf{W}^{*b}).\)

    • compute mean experience on bootstrap sample, \(\overline{\mbox{exp}}^{~*b}\).
    • run regression on bootstrap sample, obtain bootstrap sample OLS coefficients \((\widehat \beta_{0}^{*b},\widehat \beta_{1}^{*b}, \widehat \beta_{2}^{*b})\)
    • compute \[\widehat \theta^*(b) = \widehat \beta_{1}^{*b} + 2 \cdot \widehat \beta_{2}^{*b} \cdot \overline{\mbox{exp}}^{~*b}.\]

Bootstrap Replications Avg. Effect Experience

#  tab.b1 <- df[indices.b1,], bootstrap sample
mean.exp.b1 <- mean(tab.b1$experience) #  sample mean experience on bootstrap sample
reg.b1 <- lm(I(log(wage))~ experience+I(experience^2),data=tab.b1) # OLS on bootstrap sample
# form bootstrap replication of estimated mean effect experience
a <- c(0,1,2*mean.exp.b1)
theta.b1 <- as.numeric(t(a)%*%coef(reg.b1)) 

Bootstrap Replications Avg. Effect Experience

Sample Estimate, \(\widehat{\theta} =\) 0.035 \(- 2\cdot\) 0.00058 \(\cdot\) 18.14

\(=\) 0.0141

Bootstrap Replication 1, \(\widehat{\theta}^*(1) =\) 0.032 \(-2\cdot\) -0.0005 \(\cdot\) 18.05

\(=\) 0.0139

Sample Mean Exp:
18.14

\(\widehat{y}=\) 1.58 \(+\) 0.035 \(\mbox{exp}\) \(-\) 0.00058 \(\mbox{exp}^2\)

Bootstrap Mean Exp:
18.05

\(\widehat{y}=\) 1.64 \(+\) 0.032 \(\mbox{exp}\) \(-\) 0.0005 \(\mbox{exp}^2\)

\(\widehat{\theta} =\) 0.035 \(+2\cdot\) 0.0353 \(\cdot\) 18.14 \(=\) 0.0141

Bootstrap Replications Avg. Effect Experience

Sample Estimate, \(\widehat{\theta} =\) 0.035 \(- 2\cdot\) 0.00058 \(\cdot\) 18.14

\(=\) 0.0141

Bootstrap Replication 2, \(\widehat{\theta}^*(2) =\) 0.024 \(-2\cdot\) -0.00044 \(\cdot\) 17.17

\(=\) 0.0085

Sample Mean Exp:
18.14

\(\widehat{y}=\) 1.58 \(+\) 0.035 \(\mbox{exp}\) \(-\) 0.00058 \(\mbox{exp}^2\)

Bootstrap Mean Exp:
17.17

\(\widehat{y}=\) 1.72 \(+\) 0.024 \(\mbox{exp}\) \(-\) 0.00044 \(\mbox{exp}^2\)

Bootstrap Replications Avg. Effect Experience

Sample Estimate, \(\widehat{\theta} =\) 0.035 \(- 2\cdot\) 0.00058 \(\cdot\) 18.14

\(=\) 0.0141

Bootstrap Replication 3, \(\widehat{\theta}^*(3) =\) 0.023 \(-2\cdot\) -0.00034 \(\cdot\) 17.89

\(=\) 0.0107

Sample Mean Exp:
18.14

\(\widehat{y}=\) 1.58 \(+\) 0.035 \(\mbox{exp}\) \(-\) 0.00058 \(\mbox{exp}^2\)

Bootstrap Mean Exp:
17.89

\(\widehat{y}=\) 1.64 \(+\) 0.023 \(\mbox{exp}\) \(-\) 0.00034 \(\mbox{exp}^2\)

Bootstrap Replications Avg. Effect Experience

Sample Estimate, \(\widehat{\theta} =\) 0.035 \(- 2\cdot\) 0.00058 \(\cdot\) 18.14

\(=\) 0.0141

Bootstrap Replication 4, \(\widehat{\theta}^*(4) =\) 0.046 \(-2\cdot\) -0.00078 \(\cdot\) 19.19

\(=\) 0.0159

Sample Mean Exp:
18.14

\(\widehat{y}=\) 1.58 \(+\) 0.035 \(\mbox{exp}\) \(-\) 0.00058 \(\mbox{exp}^2\)

Bootstrap Mean Exp:
19.19

\(\widehat{y}=\) 1.5 \(+\) 0.046 \(\mbox{exp}\) \(-\) 0.00078 \(\mbox{exp}^2\)

Construct \(B\) Bootstrap Replications

  • We have created 4 bootstrap samples, and corresponding 4 bootstrap replications.

  • We want to create \(B\) bootstrap samples, and corresponding \(B\) bootstrap replications, for \(B\) large, for example, \(B=10,000\).

  • How to repeat steps many times?

    • we first consider approach in base R using replicate,
    • we then consider using functions from package boot.

Use replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))
  • Create user defined function and use replicate to generate \(B\) bootstrap samples, and corresponding \(B\) bootstrap replications.

Use Replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))
  • User defined function takes two arguments:
    1. original sample as data frame,
    2. indices;
  • Function returns bootstrap replication.

Use Replicate to Construct \(B\) Boot. Replications

f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
N <- nrow(df)  # Number of obs in sample
B <- 10000 # number of bootstrap replications
boot.est <-replicate(B, f.boot1(df, sample(1:N, N, replace=TRUE)))
  • Inputting user defined function, original data and random sample of indices into replicate function, which is returning \(B\) bootstrap replications.

  • Use output to compute bootstrap s.e. or CI.

Workflow using base R

flowchart LR
  A{  &nbsp  &nbsp User-Defined Fn.  &nbsp &nbsp <br>  &nbsp  &nbsp for Estimator &nbsp  &nbsp <br> &nbsp  &nbsp} --> B{   &nbsp  &nbsp  `replicate` &nbsp &nbsp}
  E{  &nbsp  &nbsp Original  &nbsp  &nbsp <br>  &nbsp &nbsp Sample  &nbsp &nbsp <br> &nbsp  &nbsp} --> B 
 B --> C{ &nbsp  &nbsp bootstrap s.e &nbsp  &nbsp <br>}
 B --> D{ &nbsp bootstrap CI &nbsp  &nbsp <br>}

Distribution of Bootstrap Estimates

Here: \[\begin{align} \theta = ~ & \beta_{1} + 2 \cdot \beta_{2} \cdot \mathbb{E}(\mbox{exp})\\ \widehat \theta = ~ & \widehat \beta_{1} + 2 \cdot \widehat \beta_{2} \cdot \overline{\mbox{exp}}\\ \widehat \theta^*(b) = ~ & \widehat \beta_{1}^{*b} + 2 \cdot \widehat \beta_{2}^{*b} \cdot \overline{\mbox{exp}}^{~*b} \end{align}\]

Key Idea:

  • Use simulated distribution of bootstrap sample statistics \(\widehat \theta^*(b)\) to approximate distribution of \(\widehat \theta\).

Distribution of Bootstrap Estimates

Graphing histogram of \((\widehat \theta^*(1), \widehat \theta^*(2), ... , \widehat \theta^*(B))\).

fig.boot <- ggplot(data.frame(boot.est), 
            aes(x=boot.est))  +
    xlab("Bootstrap Replications: Avgerage Effect Experience")+
    geom_histogram(color="darkblue", 
            fill="lightblue",binwidth=0.0005)
fig.boot 

Key Idea:

  • Use simulated distribution of bootstrap sample statistics \(\widehat \theta^*(b)\) to approximate distribution of \(\widehat \theta\).

Bootstrap Estimate Bias

Bootstrap estimate bias:

\(\bar{\theta}^* = \frac{1}{B} \sum_{b=1}^B \widehat \theta^{}(b)\)

\(\widehat{\mbox{bias}}^* = \bar{\theta}^* - \widehat \theta\)

  • Use difference between mean of bootstrap replications and estimate on original sample as bootstrap estimator of bias.

Bootstrap Estimate Bias

Bootstrap estimate bias:

mean(boot.est)-theta.hat
[1] 0.000051192
  • Use difference between mean of bootstrap replications and estimate on original sample as bootstrap estimator of bias.

    • \(\widehat{\mbox{bias}}^*=\) 0.0000512

Bootstrap Standard Errors

\[\mbox{s.e.}^{\mbox{boot}}(\widehat{\theta}_N)=\sqrt{ \frac{1}{B-1} \sum_{b=1}^B (\widehat \theta^{*}(b) - \bar{\theta}^*)^2 }\]

sd(boot.est)
[1] 0.002900513
  • Use s.d. of bootstrap replications \(\widehat \theta^*(b)= \widehat \beta_{1}^{*}(b) + 2 \cdot \widehat \beta_{2}^*(b) \cdot \overline{\mbox{exp}}^*(b)\) to form bootstrap standard errors of \(\widehat \theta = \widehat \beta_{1} + 2 \cdot \widehat \beta_{2} \cdot \overline{\mbox{exp}}\).

    • \(\mbox{s.e.}^{\mbox{boot}}(\widehat \theta)=\) 0.0029

Normal Approx. Bootstrap CI

95% normal-approximation bootstrap CI: \[C^{nb} = [\widehat{\theta}-1.96 \cdot \widehat{\mbox{se}}_B, ~\widehat{\theta}+1.96 \cdot \widehat{\mbox{se}}_B]\] Based on approximate normality.

  • Plug bootstrap s.e. into conventional 95% CI for asymptotically normal estimators.

Normal Approx. Bootstrap CI

Endpoints of Normal Bootstrap CI:

theta.hat -qnorm(.975)*sqrt(var(boot.est))
[1] 0.008424728
theta.hat +qnorm(.975)*sqrt(var(boot.est))
[1] 0.01979453
  • Plug bootstrap s.e. into conventional 95% CI for asymptotically normal estimators:
    • \(C^{nb} =\) [ 0.00842 \(,\) 0.0198 ].

Bootstrap Percentile CI

95% Bootstrap Percentile Interval: \[ C^{pc} = [q^*_{0.025}, ~q^*_{0.975}].\] where \(q^*_{0.025}\) and \(q^*_{0.975}\) are the \(0.025\) and \(0.975\) empirical quantiles of bootstrap replications.

Bootstrap Percentile CI

Bootstrap quantiles:

q <- quantile(boot.est, c(0.025, 0.975)) 
q 
      2.5%      97.5% 
0.00852697 0.01994357 
  • Bootstrap Percentile 95% CI:
    \(C^{pc} =\) [ 0.00853, 0.0199 ]

Boot. Percentile vs Norm Approx CI

  • Bootstrap Normal Approx. 95% CI:

    • [ 0.00842 , 0.0198 ].
  • Bootstrap Percentile 95% CI:

    • [ 0.00853 , 0.0199 ].
  • Why so similar?

Alternatively Can Use Package boot

Workflow using package boot.

flowchart LR
  A{  &nbsp &nbsp  &nbsp User-Defined Fn. &nbsp &nbsp &nbsp <br>&nbsp &nbsp &nbsp  for Estimator&nbsp &nbsp  &nbsp} --> B{ &nbsp &nbsp    &nbsp  `boot` &nbsp &nbsp &nbsp}
  F{ &nbsp &nbsp  &nbsp Original Sample &nbsp &nbsp  &nbsp} --> B 
 B --> C{ &nbsp &nbsp &nbsp bootstrap s.e &nbsp &nbsp  &nbsp}
  B --> D{&nbsp &nbsp  &nbsp `boot.ci` &nbsp &nbsp &nbsp}
  D --> E{ &nbsp &nbsp &nbsp bootstrap CI &nbsp &nbsp  &nbsp}

Alternatively using package boot

Function boot from package boot:

parameter description
data original sample (as vector, matrix, or data frame)
statistic function that constructs bootstrap replication on bootstrap sample, must have two arguments: (1) data, and (2) indices.
R number of bootstrap samples/replications, we have been denoting \(B\)
stype for this course, will always set to "i" for indices.

Alternatively using package boot

Function boot.ci from package boot:

parameter

description

boot.out

output from boot

conf

confidence level, by default is 0.95.

type

type of CI, includes

  • “norm” for normal approx. CI,

  • “perc” for percentile CI,

  • “bca” for bootstrap accelerated bias-corrected percentile interval.

BCa Confidence Interval

Bootstrap accelerated bias-corrected percentile interval (BCa):

  • improves upon bootstrap percentile CI by adjusting for bias and skewness in sampling distribution.

  • method not intuitive, is complicated, justification is more advanced.

  • easy to implement using boot.ci.

  • should use it when there is substantial bias/skewness in sampling distribution.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
  • user defined function takes two inputs:
    • data frame (original sample)
    • indices;
  • function returns bootstrap replication.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
  • boot taking as inputs:

    • original sample (data frame),
    • user defined function,
    • number of replications,

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = f.boot1, R = B, stype = "i")


Bootstrap Statistics :
      original        bias    std. error
t1* 0.01410963 0.00004975603 0.002904684

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
  • boot.ci taking as inputs:

    • saved output from boot,
    • choice of confidence level,
    • type of bootstrap CI.

Using package boot

library(boot)
f.boot1 <- function(d,i){ 
  mean.exp.b  <-mean(d[i,]$experience) 
  reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
  # form bootstrap estimate
  a <- c(0,1,2*mean.exp.b)
  as.numeric(t(a)%*%coef(reg.b)) }
boot.out1<- boot(data=df, statistic=f.boot1, R=B, stype="i")
boot.out1
boot.ci1 <- boot.ci(boot.out1, conf=0.95, type=c("norm","perc","bca"))
boot.ci1
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out1, conf = 0.95, type = c("norm", "perc", 
    "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   ( 0.0084,  0.0198 )   ( 0.0085,  0.0199 )   ( 0.0083,  0.0197 )  
Calculations and Intervals on Original Scale

Output from function boot

boot.out1<- boot(data=df, statistic=f.boot1, 
                 R=B, stype="i")

#boot.out1$t0 is estimate on original sample
boot.out1$t0
[1] 0.01410963
#boot.out1$t is vector of 10,000 stored 
#            bootstrap replications

ggplot(data.frame(boot.out1$t), 
            aes(x=boot.out1$t))  +
    xlab("Bootstrap Replications: 
            Avgerage Effect Experience")+
    geom_histogram(color="darkblue", 
            fill="lightblue",binwidth=0.0005)

plot(boot.out1)

summary(boot.out1$t)
       V1          
 Min.   :0.003135  
 1st Qu.:0.012088  
 Median :0.014050  
 Mean   :0.014060  
 3rd Qu.:0.015995  
 Max.   :0.024635  
summary(boot.out1)
      R original     bootBias    bootSE bootMed
1 10000  0.01411 -0.000049418 0.0028694 0.01405

Summary Avg Effect Exp

Estimated bias -0.0000494
Bootstrap s.e. 0.00287
Bootstrap Normal 95% CI [0.00837,0.0198]
Bootstrap Percentile 95% CI [0.00849,0.0199]
BCa 95% CI [0.0083,0.0197]

How to summarize bootstrap results?

What results would you report? how to interpret?

Rerun Analysis?

  • Can rerun bootstrap simulation to see if results are stable across simulations.
  • If not?
    • insufficient number of replications? increase \(B\)?
    • problem with bootstrap method in particular example?
      • sometimes bootstrap standard errors invalid/unstable, arises when bootstrap replications take extreme values, e.g., when dividing by number close to zero for some replications.

Rerun Analysis?

  • Can rerun bootstrap simulation to see if results are stable across simulations.
  • If not?
    • problem with bootstrap method in particular example?
      • sometimes bootstrap standard errors invalid/unstable, arises when bootstrap replications take extreme values, e.g., when dividing by number close to zero for some replications.
      • percentile and BCa CIs are valid under weaker assumptions and tend to be more robust/stable.

Repeated Simulations

boot.out2<- boot(data=df, statistic=f.boot1, R=B, stype="i")
plot(boot.out2)

summary(boot.out2)
      R original     bootBias    bootSE  bootMed
1 10000  0.01411 -0.000023716 0.0029324 0.014067
summary(boot.out2$t)
       V1          
 Min.   :0.002754  
 1st Qu.:0.012117  
 Median :0.014067  
 Mean   :0.014086  
 3rd Qu.:0.016044  
 Max.   :0.025135  
boot.out3<- boot(data=df, statistic=f.boot1, R=B, stype="i")
plot(boot.out3)

summary(boot.out3)
      R original     bootBias    bootSE  bootMed
1 10000  0.01411 0.0000082234 0.0029378 0.014094
summary(boot.out3$t)
       V1          
 Min.   :0.002334  
 1st Qu.:0.012077  
 Median :0.014094  
 Mean   :0.014118  
 3rd Qu.:0.016074  
 Max.   :0.026175  

Repeated Simulations

Second Simulation CIs:

boot.ci(boot.out2, conf=0.95, type=c("norm","perc","bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out2, conf = 0.95, type = c("norm", "perc", 
    "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   ( 0.0084,  0.0199 )   ( 0.0083,  0.0198 )   ( 0.0083,  0.0198 )  
Calculations and Intervals on Original Scale

Third Simulation CIs:

boot.ci(boot.out3, conf=0.95, type=c("norm","perc","bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out3, conf = 0.95, type = c("norm", "perc", 
    "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   ( 0.0083,  0.0199 )   ( 0.0085,  0.0199 )   ( 0.0084,  0.0199 )  
Calculations and Intervals on Original Scale

Summary Across Repeated Simulations

Simulation 1 Simulation 2 Simulation 3
Estimated bias -0.0000494 -0.0000237 0.0000082
Bootstrap s.e. 0.00287 0.00293 0.00294
Bootstrap Normal 95% CI [0.00837,0.0198] [0.00839,0.0199] [0.00834,0.0199]
Bootstrap Percentile 95% CI [0.00849,0.0199] [0.00828,0.0198] [0.00847,0.0199]
BCa 95% CI [0.0083,0.0197] [0.00831,0.0198] [0.00845,0.0199]

Implementing Bootstrap: Experience level that maximizes log wages.

  • Motivation: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Implementing Bootstrap: Average Effect of Experience

  • Implementing Bootstrap: Experience level that maximizes log wages.

  • Other Bootstrap Variants

Implement: Exp Level Max. Wages

  • We now implement bootstrap s.e. and CIs for experience level that maximizes log wage.

    • \(\theta = \left | \frac{\beta_{1}}{2 \cdot \beta_2} \right|,\)
    • \(\widehat \theta = \left | \frac{\widehat \beta_{1}}{2 \cdot \widehat \beta_2} \right|,\)
  • Key difference from before:

    • parameter involves ratio of coefficients,
    • estimator divides by estimated coefficient,
    • cause for worry?

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = f.boot2, R = B, stype = "i")


Bootstrap Statistics :
    original   bias    std. error
t1* 30.23942 4.237817    142.7423

Max Experience

f.boot2 <- function(d,i){ 
   reg.b <- lm(I(log(wage))~ experience+I(experience^2),data=d[i,])
   abs(coef(reg.b)[2]/(2*coef(reg.b)[3])) }
boot.out.b<- boot(data=df, statistic=f.boot2, R=B, stype="i")
boot.out.b
boot.ci.b <- boot.ci(boot.out.b, conf=0.95, type=c("norm","perc","bca"))
boot.ci.b
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out.b, conf = 0.95, type = c("norm", 
    "perc", "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   (-2013.71,  2042.13 )   (   25.29,    51.45 )   (   25.29,    51.23 )  
Calculations and Intervals on Original Scale

Output from function boot

  • Note extreme skew to the right, extreme outliers.

  • explanation?

      R original bootBias bootSE bootMed
1 10000   30.239    16.03 1034.7  30.333
       V1           
 Min.   :     2.23  
 1st Qu.:    28.09  
 Median :    30.33  
 Mean   :    46.27  
 3rd Qu.:    33.57  
 Max.   :101227.05  

Limiting X-axis:

Repeated Simulations

boot.out.b2<- boot(data=df, statistic=f.boot2, R=B, stype="i")
plot(boot.out.b2)

summary(boot.out.b2)
      R original bootBias bootSE bootMed
1 10000   30.239   4.9155 217.32  30.216
summary(boot.out.b2$t)
       V1           
 Min.   :    0.109  
 1st Qu.:   28.049  
 Median :   30.216  
 Mean   :   35.155  
 3rd Qu.:   33.457  
 Max.   :21573.185  
boot.out.b3 <- boot(data=df, statistic=f.boot2, R=B, stype="i")
plot(boot.out.b3)

summary(boot.out.b3)
      R original bootBias bootSE bootMed
1 10000   30.239   5.5227 241.72  30.257
summary(boot.out.b3$t)
       V1           
 Min.   :    2.328  
 1st Qu.:   28.015  
 Median :   30.257  
 Mean   :   35.762  
 3rd Qu.:   33.556  
 Max.   :23085.192  

Repeated Simulations

Simulation 2

boot.ci(boot.out.b2, conf=0.95, type=c("norm","perc","bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out.b2, conf = 0.95, type = c("norm", 
    "perc", "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   (-400.62,  451.27 )   (  25.14,   51.59 )   (  25.22,   53.61 )  
Calculations and Intervals on Original Scale

Simulation 3

boot.ci(boot.out.b3, conf=0.95, type=c("norm","perc","bca"))
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 10000 bootstrap replicates

CALL : 
boot.ci(boot.out = boot.out.b3, conf = 0.95, type = c("norm", 
    "perc", "bca"))

Intervals : 
Level      Normal             Percentile            BCa          
95%   (-449.04,  498.47 )   (  25.14,   52.12 )   (  25.22,   53.14 )  
Calculations and Intervals on Original Scale

Summary Across Repeated Simulations

Simulation 1 Simulation 2 Simulation 3
Estimated bias 16.03 4.915 5.523
Bootstrap s.e. 1034.7 217.32 241.72
Bootstrap Normal 95% CI [-2014,2042] [-400.6,451.3] [-449,498.5]
Bootstrap Percentile 95% CI [25.29,51.45] [25.14,51.59] [25.14,52.12]
BCa 95% CI [25.29,51.23] [25.22,53.61] [25.22,53.14]
  • what’s going on?

  • what would you report? how summarize?

Take-away lessons from application

  • We have investigated bootstrap inference for two parameters:

    • Average effect of experience: \[\beta_1 + 2 \cdot \beta_2 \cdot \mathbb{E}(\mbox{exp}),\]
    • Experience level that maximizes log wage: \[\left | \frac{\beta_1}{2 \beta_2}\right |.\]
  • How to summarize results for each parameter? What results would to rely on?

  • How are bootstrap results different for the two parameters? what do you think are the causes of the differences? what are the take-away lessons?

Other Bootstrap Variants

  • Motivation: Mincer Wage Equation

  • Delta Method vs Bootstrap

  • Bootstrap Algorithm

  • Implementing Bootstrap: Average Effect of Experience

  • Implementing Bootstrap: Experience level that maximizes log wages.

  • Other Bootstrap Variants

Paired vs Residual vs Wild Bootstrap

  • We have done nonparametric bootstrap for OLS, also called “pairs bootstrap” in regression context.

  • Other bootstrap procedures for regressions:

    • Residual bootstrap,

    • Wild bootstrap.

  • Tradeoffs?

Non-i.i.d data

  • We have done bootstrap for i.i.d. data.

  • Can modify bootstrap procedure for

    • clustered sampling
      • straight forward modification, resample clusters
    • time series
      • less straight forward (block bootstrap).

    .